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This paper shows that several known properties of the Yukawa system can be derived from the 
isomorph theory, which applies to any system that has strong correlations between its virial and 
potential-energy equilibrium fluctuations. Such “Roskilde-simple” systems have a simplified ther¬ 
modynamic phase diagram deriving from the fact that they have curves (isomorphs) along which 
structure and dynamics in reduced units are invariant to a good approximation. We show that 
the Yukawa system has strong virial potential-energy correlations and identify its isomorphs by two 
different methods. One method, the so-called direct isomorph check, identifies isomorphs numer¬ 
ically from jumps of relatively small density changes (here 10%). The second method identifies 
isomorphs analytically from the pair potential. The curves obtained by the two methods are close 
to each other; these curves are confirmed to be isomorphs by demonstrating the invariance of the 
radial distribution function, the static structure factor, the mean-square displacement as a function 
of time, and the incoherent intermediate scattering function. Since the melting line is predicted to 
be an isomorph, the theory provides a derivation of a known approximate analytical expression for 
this line in the temperature-density phase diagram. The paper’s results give the first demonstration 
that the isomorph theory can be applied to systems like dense colloidal suspensions and strongly 
coupled dusty plasmas. 


I. INTRODUCTION 

The Yukawa pair potential has been used to model 
a wide variety of different phenomena in physics [1]. 
Named after Hideki Yukawa who used this potential in 
his meson theory [2], Debye and Hiickel [3] had in fact 
used it earlier to describe the interactions between ions 
in solutions. The Yukawa pair potential, which is also 
referred to as the screened Coulomb potential, has the 
form 


n2 

u(r) = — exp(—r/A) . (1) 

r 

Here Q is the particle charge in the Gaussian unit system 
and A the so-called screening length. Although Debye- 
Hiickel theory was devised for describing the behavior of 
dilute ionic solutions, the Yukawa potential has proven 
useful also in the description of charge carriers that are 
much larger than those of the surrounding medium. This 
is the case for instance in suspensions of charge-stabilized 
colloids or in dusty plasmas, compare, e.g. the DLVO 
theory [4, 5] in which the interaction between the surfaces 
of two colloids is described by a potential of the Yukawa 
form. The colloid particles are often modeled with a 
potential that adds a hard-core repulsion, but the low- 
temperature part of the phase diagram of this “hard¬ 
core” Yukawa potential can be mapped onto the phase 
diagram of the Yukawa pair potential [6] . 
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Besides electrons and ions, dusty plasmas contain small 
solid particles that are charged. In dusty plasmas the 
distance between the solid particles is usually large com¬ 
pared to the size of the particles, and their interaction 
can therefore be modeled using the (point) Yukawa po¬ 
tential of Eq. (1). A more involved potential consisting 
of a sum of Yukawa terms is sometimes used to capture 
phenomena like plasma production and loss balance [7]. 
The physics of dusty plasmas is not only of interest for 
industrial applications where such plasmas are formed, 
but also in astrophysics for the understanding of stellar 
materials and planet formation [8, 9]. 

In 1986 it was found that colloids in suspension can 
form crystal lattices [10] and it was predicted that this 
should also be possible in dusty plasmas [11]. This cre¬ 
ated a renewed interest in the phase diagram of the 
Yukawa system [12-14]. It was later found in experi¬ 
ments that plasma crystals indeed do exist [15, 16]. If p 
is the number density of particles, in the field of dusty 
plasmas the phase diagram is usually presented in terms 
of the following two dimensionless parameters [17, 18]: 
the screening parameter 



and the coupling parameter 


keT 


( 2 ) 

( 3 ) 


Physically, the screening parameter is much larger than 
unity whenever the screening length is much smaller than 
the average interparticle distance - this characterizes the 
part of the thermodynamic phase diagram where the ex¬ 
ponential damping term dominates over the Coulomb 
term. In the other limit, k <C 1, the exponential damping 
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that some model liquids have strongly correlated fluctu¬ 
ations of their energy and pressure [34]. More precisely, 
the correlations are between the configurational parts of 
these quantities. If one for a system of N particles splits 
energy and pressure into ideal-gas terms, which depend 
only on the particle momenta, and configurational terms 
depending only on the particle positions as follows 

E = K{pi,...,pN) + U{ri,...,rN), (4) 

pV = NksTipi,.. .,Pn) + W{ri,...,rN), (5) 


FIG. 1. Scatter plot (red symbols) demonstrating strong cor¬ 
relations between the instantaneous values of the virial W and 
the potential energy U during an equilibrium NVT simulation 
of the Yukawa system at p = 3 x 10~® and T = 5.3608 x 10“® 
(k = 6.934 and F = 2690). The dashed line is the standard 
linear regression with slope 7 . 


term plays little role and the system behaves as a single¬ 
component repulsive Coulomb system, the so-called one- 
component plasma (OCP), with an ill-defined thermo¬ 
dynamic limit due to the infinite screening length. The 
case r 1 corresponds to the potential energy from 
the individual pair interactions being much larger than 
the thermal energy; this favors crystallization, depend¬ 
ing on the value of k, see below. In the opposite limit, 
r <g; 1, the system approaches an ideal gas. Care must 
be taken when comparing k or F values in the litera¬ 
ture, because they often include in their definition fur¬ 
ther multiplicative constants. For instance, in the above 
definitions of the screening and coupling parameters the 
term is often replaced by the Wigner-Seitz radius 

a = [3/(47rp)]i/3. 

The Yukawa system’s phase diagram is well under¬ 
stood with its two solid phases and a fluid phase [11- 
14]. At high K the Yukawa fluid crystallizes into an 
face-centered cubic (FCC) lattice, while at low k it crys¬ 
tallizes into a body-centered cubic (BCC) lattice. The 
triple point separating fluid, FCC, and BCC phases is 
located at k = 6.90,F = 3.47 x 10^ [19]. Vaulina et 
al. [20, 21] derived an expression for the melting line from 
the pair potential using the Lindemann melting criterion 
[22, 23]. This expression is confirmed below, where it 
is derived from the isomorph theory. This theory ap¬ 
plies to the class of systems (dense liquids and classical 
crystals) termed Roskilde-simple - or just Roskilde {R) 
- systems [24, 25]. As documented below, the Yukawa 
system exhibits the strong correlations between equilib¬ 
rium virial and potential-energy fluctuations required for 
a system to belong this class. 

Over the years, investigations of the dynamics of the 
Yukawa system have mostly focused on the self-diffusion 
coefficient and the viscosity [12, 17, 26-33]. In the present 
publication we derive some of the previously reported 
scaling properties of the Yukawa system’s transport co¬ 
efficients from the general isomorph theory. 

The isomorph theory initiated from the observation 


strong correlations between the potential energy U and 
virial W constant-volume thermal-equilibrium fluctua¬ 
tions are found for a number of systems [35] . An illustra¬ 
tion of such strong virial potential-energy correlations is 
shown in Fig. 1 that plots the instantaneous values of W 
versus those of U during an NVT equilibrium simulation 
of the Yukawa system. The correlations are quantified by 
the standard Pearson correlation coefficient [36] 
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As Fig. 1 shows, R is high for the Yukawa system. A 
system is generally considered to have strong C/, W cor¬ 
relation whenever R > 0.9 [36], a convenient but also 
somewhat arbitrary criterion. Many model systems have 
been shown to exhibit strong U, W correlations [36-38] , 
including the Lennard-Jones (LJ) system and other sim¬ 
ple liquids, molecular liquid models [39, 40], liquids under 
shear [41] , and crystals [42] . These systems were initially 
called “strongly correlating”, but that repeatedly led to 
confusion with strongly correlated quantum systems, and 
they are now instead referred to as Roskilde-simple or 
just Roskilde (R) systems [43-50]. The isomorph theory 
has also been applied to nano-confined liquids [51, 52] 
and, e.g., to derive density-scaling invariants for zero- 
temperature plastic flow properties of glasses [53]. 

In 2009 it was found that R liquids have “isomorphic” 
curves in their thermodynamic phase diagram [54]. Iso- 
morphs are curves along which all properties derived from 
structure or dynamics are invariant in properly reduced 
units, making the phase diagram effectively one dimen¬ 
sional with respect to many properties [54] . An example 
of this is shown in Fig. 2, where the Yukawa fluid is shown 
to have - to a very good approximation - identical radial 
distribution functions (<?(?’)) and mean square displace¬ 
ments at two state points that are isomorphic to each 
other. 

Having its origins in the field of organic, glass-forming 
liquids, the isomorph theory has been applied to explain 
certain empirical scalings in that field. These include the 
so-called density scaling, according to which the dynam¬ 
ics is a function of p'^/T [55]. This is an approximation 
to the isomorph theory that applies when fairly small 
changes in density are considered, as in most experiments 
[54, 56]. Because both the excess entropy, Sgx = S—Sideai 
{Sideai being the entropy of the ideal gas at same tem¬ 
perature and density), and the reduced-unit dynamics 
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FIG. 2. The radial distribution function (left) and the mean 
square displacement (right) in reduced units for four differ¬ 
ent state points, two of which are isomorphic to each other. 
The initial state point (black line) has pi = 3 x 10~®,ri = 
5.361 X 10“®, while the isomorphic state point (dashed, red 
line) has p 2 = 3.3 x 10 ~®,r2 = 6.553 x 10~®. The method 
used to identify the state point (p 2 , Tz) as being isomorphic to 
state point (pi,ri) is given below in Sec. IV A. The two iso¬ 
morphic state points are seen to have identical structure and 
dynamics to a very good approximation. The effects of ei¬ 
ther changing density or temperature (green and blue dashed 
lines) are shown for comparison. The insets show enlarged 
views of the same data. 


are invariant along isomorphs, the isomorph theory also 
predicts that R liquids obey Rosenfeld’s excess entropy 
scaling, which states that the reduced transport coeffi¬ 
cients are functions of Sex [57]. Not predicted by the 
isomorph theory, but nevertheless found to apply to R 
liquids [58], is the Rosenfeld-Tarazona relation [59] that 
gives expressions for the constant-density temperature 
dependence of the isochoric specific heat (cy oc 
For a more comprehensive overview of the isomorph the¬ 
ory the reader is referred to a recently published Feature 
article [25]. 

The fact that the Yukawa system has strong UW cor¬ 
relations shows that it belongs to the class of R liquids. 
This is consistent with the long-known facts that the 
Yukawa system [28, 60] and the OCR [61] obey excess 
entropy scaling as well as the above-mentioned Rosenfeld- 
Tarazona relation for the specific heat’s temperature de¬ 
pendence [59] . The implication is that the isomorph the¬ 
ory applies also to systems that are not conventional liq¬ 
uids, such as colloidal suspensions and dusty plasmas. In 
fact, the findings presented below confirm a recently de¬ 
rived theory of quasiuniversality according to which the 
Yukawa system is in the so-called EXP quasiuniversal¬ 
ity class of simple, monatomic, pair-potential systems, a 
class that includes also, e.g., the inverse-power-law (IPL) 
and Lennard-Jones systems [62]. 

In the following we first give a short overview of the 
isomorph theory and the computer simulation procedures 
used in the paper (Secs. II and III). We proceed to show 
how to construct isomorphs using two different methods 
in Sec. IV. Sec. V then shows that the isomorphs gener¬ 
ated by these methods are very close to one another and 
confirms the predicted isomorph invariance of structure 


and dynamics. Finally, Sec. VI gives a brief discussion. 


II. ISOMORPHS IN ROSKILDE-SIMPLE 
SYSTEMS: A BRIEF REVIEW 

The isomorph theory uses so-called reduced units in 
which quantities are made dimensionless via units based 
on the thermodynamic quantities density and tempera¬ 
ture, not the length and energy of the microscopic poten¬ 
tial that are often used in reporting simulation results. 
Thus one uses the unit Iq = for length, eo = ksT 

for energy, and to = y/ksT/m for time with m 

being the particle mass (for Brownian dynamics a differ¬ 
ent time unit is used) [54]. For instance, the collective 
position of the system’s particles, R = (ri,..., r^r), is 
expressed in reduced units as 

R = (7) 

where the tilde here and henceforth denotes a reduced 
quantity. 

The isomorph theory is conveniently summarized in 
the expression [25] 

17(R)-h(p)$(R)+5(p). (8) 

Here $ is a dimensionless, state-point independent func¬ 
tion of the reduced configurations R; the term $(R) 
controls the structure and dynamics in reduced units, 
but notably it contains no length or energy scales. The 
functions h{p) and g{p) both have units of energy. The 
physics of Eq. (8) is that upon changing the density of 
a system, the potential-energy surface to a good approx¬ 
imation merely undergoes a linear, affine rescaling. We 
proceed to show how Eq. (8) is used to derive the most 
important properties of R systems, their strong virial en¬ 
ergy correlations and the existence of isomorphs. Before 
doing so we note that the isomorph theory was recently 
generalized by defining a Roskilde-simple system by the 
property that the order of potential energies is main¬ 
tained for uniform scaling of configurations, i.e., by the 
condition C/(Ra) < C/(Rb) => U{XRa) < 17(AR{,) [63]. 
This leads to slightly modified predictions, but overall 
the new isomorph theory is close to the original, which 
is the one used below. 

The microscopic virial of Eq. (5) is defined by [64] 

W(R) = -iR-VC7(R). (9) 

It is easy to show that kF(R) characterizes the change 
of potential energy for a uniform scaling of space (i.e., 
leaving R intact) as follows [65] 
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Combining this with Eq. ( 8 ) and eliminating $(R) one 
finds that 

VE(R) ^ 7(p)t/(R) + <^(P) • (11) 

Here 


the isochoric specific heat, as well as the reduced-unit 
dynamics and structure. Not all quantities are invariant 
even in reduced units, for instance the free energy and 
its volume derivatives like pressure or compressibility. 

Because the temperature is proportional to h{p) along 
an isomorph with the same proportionality constant, iso- 
morphs are described by 


,n, 

and 4>{p) = dg{p)/ p — g{p)^{p). From Eq. (11) it 

follows that at constant density the fluctuations in W 
and U are correlated with linear-regression slope 7 given 
by 


—^ = Const. (16) 

Since moreover the excess entropy per particle Sex is con¬ 
stant on an isomorph (because the Boltzmann probabili¬ 
ties of scaled configurations are), the temperature of an R 
system, i.e., one with isomorphs, separates as follows [ 66 ] 


_ (AlEAC/) 

((At/) 2 ) 


(13) 


ksT = f{sex)h{p) . (17) 


As shown in Fig. 1(a) the Yukawa system indeed has 
strong WU correlations (with slope 7 = 2.14 at the state 
point studied here). 

We proceed to show that Eq. ( 8 ) implies the exis¬ 
tence of isomorphs. First these curves in the thermo¬ 
dynamic phase diagram need to be defined. Consider 
two densities, pi and p 2 , and two configuration Ri and 
R 2 at these densities, respectively, with the same re¬ 
duced coordinates, i.e., pJ^^Ri = = R- Ap¬ 

plying Eq. ( 8 ) to 17(Ri) and C/(R 2 ) and eliminating the 
common factor $(R) we find [C/(Ri) — g{pi)] /h{pi) = 
[C/(R 2 ) — 3 (^ 2 )] /h{p 2 )- Define now two temperatures by 
ksTi = Kh{pi) and kBT 2 = Kh{p 2 ) with the same pro¬ 
portionality constant K. This gives 


B(Ri) ^ £/(R2) , „ 

fcsTi kBT2 


(14) 


where the constant B 12 depends only the two state 
points, not on the configurations. This can be rewrit¬ 
ten as 


exp 


u(Ri) \ 

kBTi ) 


C \2 exp 


( U(R2) \ 

V kBT2 ) ■ 


(15) 


This equation implies identical canonical probabilities of 
configurations with the same reduced coordinates. When 
this is obeyed to a good approximation for most of the 
physically relevant configurations, the two state points 
(pi,Ti) and (^ 2 , 22 ) are by definition isomorphic to one 
another [54]. This defines a mathematical equivalence 
relation in the thermodynamic phase diagram, and an 
isomorph is then defined as a curve along which all pairs 
of state points are isomorphic. 

From the isomorph definition Eq. (15) it follows that 
many properties of the system are invariant between iso¬ 
morphic state points [54]. Isomorph invariants include 
thermodynamic quantities such as the excess entropy and 


It can be shown that this separation property is mathe¬ 
matically equivalent to the thermal average of Eq. (11), 
which is simply the configurational Griineisen equation 
of state W = 'y{p)U + (j){p) [ 66 ]. 

The isomorph theory is approximate for all realistic 
systems. In fact, Eqs.( 8 ) and (15) are only exact for sys¬ 
tems that have an Euler-homogeneous potential-energy 
function. The most important example of such systems 
are those with an inverse-power-law (IPL) pair potential 
{v{r) oc r“”). IPL systems are easily shown to be char¬ 
acterized by h{p) oc p'*', and it follows from Eq. (9) that 
7 is related to the exponent n of the IPL pair potential 

by 



The applicability of the isomorph theory for simple 
atomic systems may be understood physically from the 
fact that, in the region of the first peak of the radial dis¬ 
tribution function, the interatomic potential is well fitted 
by an IPL term plus a linear term, the so-called extended 
IPL (eIPL) pair potential [67]. This provides an intu¬ 
itive explanation why many non-IPL liquids obey the 
isomorph theory. Notably, the linear term contributes 
little to the fluctuations in the energy or the virial. The 
reason for this is that when a particle moves, interparti¬ 
cle distance decrease on one side while increasing on the 
other side. The sum of the interparticle distances stays 
approximately the same. [67] 


III. SIMULATION PROCEDURE 

Our simulations used the “Roskilde University Molec¬ 
ular Dynamics” (RUMD) code, which is optimized for 
GPU computing [ 68 ]. All simulations were performed in 
the NVT ensemble with a standard Nose-Hoover thermo¬ 
stat. A state point is characterized by two parameters: 
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N 

^ cut 

At 


p < 0.01 

2048 

4.3 

0.0010 

5 < K 

0.01 ^ p < 0.5 

2048 

5.2 

0.0025 

1.5 < K < 5 

0.5 <p 

8192 

10.0 

0.0025 

re < 1.5 


TABLE I. The parameters used in the simulations. Depend¬ 
ing on the density p the number of particles in the system 
N was changed to allow for a larger reduced cutoff radius 
fcut- The reduced time step At was decreased at the lowest 
densities. 


the density p, reported below in the unit system defined 
by the screening length A of Eq. (1), and the tempera¬ 
ture T reported in the unit system defined by the unit 
Q'^/{Xkg). Equivalently, a state point may be charac¬ 
terized by the two dimensionless numbers the screening 
parameter k of Eq. (2) and the coupling parameter E 
of Eq. (3). Eor maximum clarity the simulation results 
are presented both in terms of p and T, and in terms 
of K and r. In practice the simulations were performed 
directly in reduced units, i.e., density and temperature 
were set equal to unity, changing instead the length and 
energy parameters of the potential in order to investi¬ 
gate different state points in the thermodynamic phase 
diagram. 

In reduced units the integration time step At was de¬ 
creased at low screening lengths (/c » 1) because the 
potential here becomes very steep. We used a shifted- 
force cutoff [69] for the potential, with a (reduced) cutoff 
fcut that varied with the state point because longer cut¬ 
offs are needed at low k values. Wherever necessary we 
increased the number of particles N to allow for larger 
cutoffs. It should be noted that this cutoff method does 
not yield accurate values of the potential energy, espe¬ 
cially at high densities close to the OCP limit. The accu¬ 
rate calculation of the potential energy of the OCP hase 
been discussed elsewhere [70-72], and is beyond the scope 
of this investigation. Table I summarizes the simulation 
parameters used. 


IV. TWO METHODS FOR IDENTIFYING THE 
ISOMORPHS 


This section details two methods for mapping out an 
isomorph in the thermodynamic phase diagram. The 
first one is numerical; it is accurate for small density 
variations, but not well suited for studies involving large 
density variation. Here a recently proposed approximate 
analytical method is more handy. This method estimates 
the isomorph from the pair potential by using Eq. (16) in 
conjunction with an analytical expression for h{p). This 
has been shown to work well for the Lennard-Jones sys¬ 
tem and two systems with pair potentials that are sums 
of two, respectively three IPL terms [73]. 



pi = 0.003 
p2 = 0.0033 


R = 0.99991 
TjTi = 1.2224 


FIG. 3. The Yukawa system simulated at pi = 3 x 10”^ 
(k = 6.93) and Ti = 5.361 x 10“® (T = 2690). The potential 
energies U\ = [/(Ri) of configurations were plotted versus 
the potential energies U 2 = [/(R2) of the same configurations 
scaled to a 10% higher density denoted p 2 (green symbols). 
The two energies are highly correlated, indicating that the 
Yukawa system obeys the isomorph definition (Eq. (14)). By 
simple linear regression (dashed line) we find the slope T^jTx 
(see Eq. (19)), which means that for the state point at pi to 
be isomorphic to the one at p 2 , the temperature of the latter 
should be T 2 = 1.2224 x Ti = 6.553 x 10“® (T = 2272). 


A. The direct isomorph check: A numerical 
method 

Equation (14) can be rewritten as 

[/(R2) = ^C/(Ri) + Di2. (19) 

4l 

This shows that the potential energies of configurations 
at density pi and the same ones scaled to density p 2 are 
predicted to be linearly related. If the configurations 
are taken from an equilibrium simulation at state point 
(pi,Ti), the slope is the ratio T 2 /T 1 , where T 2 is the 
temperature for which state point (p 2 ,T 2 ) is isomorphic 
to state point (pi,ri). This can be used to find state 
points that are isomorphic to each other. This procedure 
is termed the direct isomorph check because it is based 
directly on the isomorph definition Eq. (15) [54]. 

The direct-isomorph-check procedure is illustrated in 
Eig. 3. Here a number of configurations Ri with poten¬ 
tial energy denoted by C/i were sampled from an equi¬ 
librium simulation of the Yukawa fluid at density pi. 
Each configuration was then scaled to the higher den¬ 
sity P 2 = 1-1 X Pi, at which the new potential energy 
U 2 was calculated. The figure shows that the potential 
energies Ui and U 2 are highly correlated. The standard 
linear-regression slope is 1.2224. According to Eq. (19) 
this number is the ratio T 2 IT 1 , which allows for determin¬ 
ing T 2 such that the state point (p 2 ,T' 2 ) is isomorphic to 
(pi, Ti). We see that along the isomorph through (pi, Ti) 
a density increase of 10% implies a 22.2% increase in tem¬ 
perature. 

The structure and dynamics of the two isomorphic 
state points are shown to be the same in Eig. 2, indicat¬ 
ing that the direct isomorph check works well for density 
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FIG. 4. Estimates of the effective IPL exponent for 

the Yukawa potential using Eq. (20). 


changes of 10%. However, it should be noted that the 
change in density cannot be huge, because the isomorph 
theory is only approximate (except for IPL systems), 
which implies that direct-isomorph-check plots give rel¬ 
atively poor correlations for large density jumps. To 
avoid this problem, we used in this work always a density 
change of merely 10 % for direct isomorph checks, corre¬ 
sponding to changing n by less than 4%. In the simula¬ 
tions presented in section V we used the direct isomorph 
check to create an isomorph by doing a simulation at the 
initial state point (p,T) = (10“^,3 x 10“®), scaling con¬ 
figurations to a new density, and finding the temperature 
of the new isomorphic state point. This was repeated at 
the new state point to obtain a third one, etc. In this way 
we obtained altogether a set of 32 prospective isomor¬ 
phic state points in the range 3 x 10“^ ^ p ^ 3.6 x 10“® 
(6.5 ^ K ^ 14.9). These state points are referred to here 
as “prospective” because they will be compared to other 
sets of prospective isomorphic state points generated in 
different ways. 


B. Predicting the isomorph analytically from the 
pair potential 

In this section we aim to construct an isomorph by 
obtaining an expression for the function h[p) of Eq. ( 8 ). 
For pair potentials that are a sum of IPLs, h{p) can be 
determined from a single simulation at a reference state 
point [56, 66 ]. For other potentials such as the Yukawa 
this is not possible. Nevertheless, Bphling et al. [73] have 
recently shown that h{p) can be estimated from the po¬ 
tential. We briefly review these findings before applying 
this method to the Yukawa potential. 

As mentioned earlier, for an IPL pair potential oc r“" 
it is known that 7 = n/3 (Eq. (18)). For other potentials 
an effective r-dependent IPL exponent can be estimated 
using ratios of derivatives of the potential [67] as follows: 

n(P\r)^-r - p, (20) 

where denotes the pth derivative of the potential. 
For an IPL pair potential, n^P'^{r) is constant and gives 


the correct exponent for all p. For other potentials, the 
effective exponent depends on both p and r, meaning that 
the “softness” of the particles depends on the separation 
between the particles. We show this in Fig. 4 for the 
Yukawa potential. It is known that in the OCP limit the 
potential reduces to a Coulomb interaction, for which 
n = 1. This corresponds to small interparticle distances 
and we see indeed that in the limit of r —> 0 one hnds 
77 ,(p) ( 7 -) X for every p. 

Recall that potentials that obey the isomorph theory 
have been found to be well fitted by an extended IPL 
potential, i.e., an IPL potential plus a linear term [67]. 
An obvious choice for p is thus 2, since this would ignore 
the linear contribution to the extended IPL potential: 


'nP\r) = 



( 21 ) 


The next question is at which distance to evaluate 
Since the physics of R liquids has been shown 
to be governed by the interactions in the first coordina¬ 
tion shell [24], the nearest neighbor distance is an ob¬ 
vious choice. Since any distance scales with density as 
r oc we can write the nearest-neighbor distance as 

Kp~^/^, where A is a number close to unity. By combin¬ 
ing Eq. (21) and Eq. (18) one finds [73] 

( 22 ) 



As noted by Bphling et al. [73], A is not expected to be 
the same for all state points. However, because structure 
in reduced units is predicted to be invariant along an iso¬ 
morph, A must also be an isomorph invariant. By com¬ 
paring with simulation results of different pair potentials, 
they indeed find that 7 (p) does not predict the density¬ 
scaling exponent equally well for different isomorphs if 
A is isomorph independent, although the prediction is 
qualitatively correct. Instead, an isomorph dependent 
(and thus excess entropy dependent) nearest-neighbor 
distance A{sex)p~^^^ was able to estimate the isomorph 
more precisely, giving [73] 


lip, Sex) 


(r) 


r=A(Sea,)p-l/3 


(23) 


It is possible to rewrite Eq. (21) as rS^\r) = 
dln[r^t>"(r)]/dlnr. This can be combined with Eq. (12) 
to find [73] 

hip, Sex) = , ( 24 ) 

where A is an arbitrary constant. This can now be used 
to find an isomorph in the phase diagram using the facts 
that h(p, Sex)/T = Const. (Eq. (16)) and Sex is an iso¬ 
morph invariant. 

Applying this method to the Yukawa potential, we find 
using Eq. (24) after straightforward calculations 


hip, Sex 




Ap“^/®-I-2-I-2A“ ■ 

(25) 
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If we express h{p, Sex) in terms of the Yukawa parameters 
K and r, we find 


T. _ T- 

^ ^ (Ak)2 + 2Ak + 2 


(26) 


where Fq = A/(2A) is the value of F in the OCP limit 
(k = 0). Vaulina et al. [20] found that a curve of this 
shape with A = 1 and Fq = 106.6 gives a good descrip¬ 
tion of the melting line as found by Hamaguchi et al. [19], 
and also of the melting line of dissipative Yukawa sys¬ 
tems [21]. Their findings are fully consistent with the iso- 
morph theory and now placed in a more general setting: 
The melting line is predicted to be an isomorph [54] since 
if an isomorph were to cross the melting line this would 
mean Eq. (15) should be obeyed for a pair of liquid and 
solid state points. This cannot be the case, because the 
relative Boltzmann probabilities for liquid and solid con¬ 
figurations are clearly different at these two state points. 
Vaulina et al. used Lindemann’s melting criterion which 
states that melting happens when the root mean square 
displacement is 10 % of the crystal’s nearest neighbor dis¬ 
tance [23]. Since both structure and dynamics are invari¬ 
ant on the isomorph, they should also be so at the melt¬ 
ing line. Lindemann’s melting criterion is thus consistent 
with the isomorph theory [54]. It should also be noted 
that Yazdi et al. [18] recently showed that Eq. (26) with 
Fq = 368 and A = 1 gives a good description of the ideal 
glass transition of mode-coupling theory; apparently this 
line is an isomorph which makes sense since isomorphs 
are lines of identical physics (in reduced units). A con¬ 
sequence of the melting line being an isomorph is also 
that the temperature T on an isomorph scaled by the 
melting temperature is an isomorph invariant if iso¬ 
morphs have the same shape (i.e., their shape is given 
by h{p) instead of h{p,Sex))- Thus so-called melting 
temperature scaling methods that have found that dy¬ 
namic and/or thermodynamic properties of the Yukawa 
system [28, 29, 75, 76] and the OCP [28, 29, 77] are func¬ 
tions only of T/Tm, are saying that these properties are 
isomorph invariants. 

In addition to the direct isomorph check described in 
section IV A we used Eq. (25) to obtain two sets of iso¬ 
morphic state points. For this it is necessary to know 
the relevant reduced interparticle distance A. Because 
Vaulina et al. [20] found Eq. (26) with A = 1 to be a 
good description of the melting line and the melting line 
is an isomorph, we test if Eq. (25) with A = 1 is the 
correct description of an isomorph. The prospective iso¬ 
morph with A = 1 that we have tested is characterized 
by A = 6.347 x lO'^ (Fq = 78.78). 

We also tried to calculate a more accurate value of A. 
Previously, Bphling et al. [73] used the most probable 
nearest neighbor distance for A, which they determined 
from the position of the first peak in r^g(r), with g(r) 
being the radial distribution function. They found for 
the potentials they tested at different state points that 
0.975 < A ^ 1.065. We below use the different method to 
determine A from the U, W fluctuations published earlier 


P 



FIG. 5. Prospective isomorphs in the (a) p, T phase diagram 
and the (b) k, F phase diagram generated in the three differ¬ 
ent ways described in Sec. IV. Comparing the two estimates 
of h{p) from the pair potential, there is a visible difference 
in the slope at low density (high k). The isomorphs are ap¬ 
proximately parallel to the expressions of the melting line [20] 
and the ideal MCT glass transition [18] (dashed lines), show¬ 
ing that these lines are both isomorphs. This is what one 
expects from the fact that the physics is invariant along the 
isomorphs. Sets of isothermal and iso-F state points have also 
been included in the figure because these are used later in this 
section. 


by Bailey et al. [74]. First, we derive the expression for 
lip, Sex) = d\nh{p,Sex)/d\np (Eq. (12)) for the Yukawa 
potential using Eq. (25): 


7(p, Sex) 3^2^173 + 0Ap2/3 + p + 3 ' 

We then did a simulation at an initial state point {p = 
0.005 and T = 0.00015). At this state point 7 was found 
from the fluctuations to be 1.78 using Eq. (13). From 
Eq. (27) it then follows that A = 1.03 is the relevant re¬ 
duced nearest-neighbor distance at that state point. The 
third prospective isomorph is thus identified by Eq. (25) 
with A = 1.03 and A = 7.414 x lO-^ (Fq = 69.46) to 
obtain a set of state points with the same densities as 
the points generated for the A = 1 prospective isomorph. 


V. COMPARING THE THREE PROSPECTIVE 
ISOMORPHS 

This section compares the three sets of prospective iso¬ 
morphic state points generated as described in the previ¬ 
ous section. The three prospective isomorphs are shown 
in Fig. 5 plotted both in the p, T and the n, F plane. The 
isomorphs are parallel to the expression for the melting 
line [ 20 ] and the ideal glass transition line from mode 
coupling theory [18]. The three prospective isomorphs 








1.00 


oi 


?- 




FIG. 7. Relative difference between the estimated 7 (p) 
(Eq. (27)) and 7 calculated from the fluctuations via Eq. (13) 
(denoted in the figure by 7 /). The prediction 7 (p)a=i .03 is 
generally in best agreement with the values from the equilib¬ 
rium fluctuations. 



FIG. 6 . The virial potential-energy correlation coefficient R 
(top), and the so-called density-scaling coefficient 7 (bottom) 
plotted versus density (left) and the screening parameter k 
(right). The Yukawa fluid is Roskilde-simple at all simulated 
state points. We find that 7 varies significantly, reflecting a 
dramatic change in steepness of the Yukawa potential with 
density. 


are slightly different, but overall close to one another. It 
is therefore not surprising that they have similar invari¬ 
ance properties. The remainder of this section focuses on 
illuminating the minor differences. 

Figures 6 (a) and (b) verify that the strong correlations 
between the virial and the potential-energy NVT equi¬ 
librium fluctuations are present over the entire range of 
densities studied. Whereas models of conventional liq¬ 
uids are usually not strongly correlating at low densities 
and temperatures due to the large contribution of the 
attractive term in the potential, the Yukawa system is 
in fact strongly correlating with correlation coefficient 
R > 0.99 (Eq. ( 6 )) at all tested state points. The val¬ 
ues R > 0.99 are very high, especially when compared to 
those of other models that, as mentioned, are considered 
Roskilde liquids if they obey R > 0.90 [36, 37, 39, 67]. 
Recently, the ten-bead rigid-bond flexible Lennard-Jones 
chain was even demonstrated to have excellent isomorphs 
despite having R ^ 0.86 [40]. 

Figures 6 (c) and (d) show the density-scaling coeffi¬ 
cient 7 calculated using Eq. (13). In the OOP limit we 
hnd that 7 goes to 1/3 as expected for an IPL with n = 1. 
At lower densities the behavior of 7 is inhuenced by the 
exponential term, for which the estimated 7 (/o) has also 
been shown for comparison. Due to the increase in the 
steepness of the potential, 7 increases from 1/3 to 5 in 
the density range shown. The lowest correlation coeffi¬ 
cient is found in the crossover region of densities marking 
the region where the effects of both the exponential and 
Coulomb terms are important. 

There is a slight difference between the values of 7 
calculated from the huctuations and those predicted for 


FIG. 8 . (a) The reduced mean-square displacement and (b) 
the self intermediate scattering function plotted as functions 
of the reduced time on the isomorph generated using the di¬ 
rect isomorph check. The dynamics collapse almost perfectly. 
The q vector is kept constant in reduced units as q = p^^^q 
with q = 7. 


"yip, Sex) from Eq. (27) with A = 1 (black dashed line). 
To investigate this further we plot in Fig. 7 the relative 
difference between the two estimated functions 7 (p, Sex) 
for A = 1 and A = 1.03 and the values of 7 calculated 
from the huctuations. Both predictions are too low at 
high density (small k). Unsurprisingly, the prediction 
with A = 1.03, a number that was identihed using the 
value of 7 from the energy and virial huctuations at a 
reference state point with p = 0.005 (k = 5.85), is more 
accurate at low densities (large k). 

We proceed to check to which degree the three sets of 
prospective isomorphic state points exhibit the invariance 
of dynamics and structure predicted for isomorphs. We 
test the invariance of the dynamics in Fig. 8 for the iso¬ 
morph generated by the direct-isomorph-check method 
and in Fig. 9 for the two isomorphs generated from the 
estimates of hip). The hgures show the mean-square dis¬ 
placements and the intermediate scattering functions in 
reduced units. For all three sets of prospective isomor¬ 
phic state points both measures of the dynamics collapse 
nicely onto a single curve. 

There are only minor differences in how invariant the 
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FIG. 9. (a) and (c) The reduced mean-square displacements 
(left); (b) and (d) intermediate scattering functions (right) 
along the two isomorphs obtained from Eq. (25) for A = 1 
(top) and A = 1.03 (bottom). Both sets of state points show 
almost invariant dynamics, although for A = 1.03, the col¬ 
lapse is slightly worse than for A = 1. The reduced scattering 
vector was again q = 7. 
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FIG. 10. The reduced diffusion coefficient along the three 
prospective isomorphs, calculated from linear regression of 
the part of the mean-square displacement above unity. The 
diffusion coefficients of the state points obtained with the di¬ 
rect isomorph check are constant in the tested range of den¬ 
sities. The two estimates from Eq. (26) have been tested in 
much larger range, and show a larger deviation from isomorph 
invariance, especially the estimate with A = 1.03 at high den¬ 
sity. Overall, the dynamics on either of the three prospective 
isomorphs are much more invariant, however, than on the 
curves of constant temperature or constant F. 


dynamics are on the three prospective isomorphs. To 
amplify these differences we calculated the reduced dif¬ 
fusion coefficient from the mean-square displacement at 
all investigated state points. The results are shown in 
Fig. 10, and compared with the diffusion coefficient along 
an isotherm and a curve of constant T. In view of the 
large density range simulated there is little variation in 
the diffusion coefficient along the isomorph. The iso¬ 
morph obtained with the direct-isomorph-check method 


FIG. 11. The reduced radial distribution function g{r) at iso¬ 
morphic state points obtained with the direct isomorph check 
(a). Although the change in density is an entire decade, the 
reduced-unit structure does not change significantly along the 
isomorph. The invariance is especially striking when com¬ 
pared to the reduced radial distribution functions of isother¬ 
mal state points (b), especially considering the smaller density 
range for the isotherm. 


only covers part of the phase diagram (because of the 
time-consuming nature of obtaining these state points). 
The reduced diffusion coefficients are virtually constant 
over the range simulated by this method. The results 
for the isomorph estimated from the pair potential with 
A = 1.03 collapse with those of the isomorph from the 
direct isomorph check. Interestingly, the range where the 
diffusion coefficients are almost invariant coincides with 
the range where agreement between the estimated and 
the fluctuation 7 is best. We also note that although 
A = 1 leads to a worse prediction for 7, it gives a more 
invariant diffusion coefficient when the whole phase dia¬ 
gram is considered. We attribute this to some cancella¬ 
tion of errors, rather than reflecting that A = 1 gives a 
more precise value of the relevant interparticle distance. 
Note also that the invariance of the dynamics, and thus 
the isomorphs seems to continue all the way to the OCP 
limit. As mentioned earlier, also the melting line is an 
isomorph, so our results indicate that also the melting 
line continues to the OCP limit, indicating the existence 
of a phase transition in the OCP [78], at least from a 
dynamical point of view. 

We proceed now to testing the isomorph invariance of 
structure as quantified via the radial distribution func¬ 
tion g{r) and the static structure factor S{q) (both in 
reduced units). Results for the prospective isomorphic 
state points obtained using the direct isomorph check 
are shown in Fig. II; results for state points of the two 
estimates of h{p) with A = I and A = 1.03 are plotted in 
Fig. 12. Overall, the structure is invariant on the three 
prospective isomorphs when compared to the change in 
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FIG. 12. The radial distribution functions (left) and static 
structure factors (right) for state points obtained from 
Eq. (25) for two values of A. The structure is not completely 
invariant on the isomorph; both the height and the position 
of the first peak vary somewhat. 

g{f) along the isotherm. This is especially impressive 
considering the density on the isotherm changes by a fac¬ 
tor of 4.4, while on the isomorphs it changes by a factor 
of 10 (Fig. 11(a)) and 10"^ (Fig. 11). 

Whereas the 10% density change applied in Fig. 2 re¬ 
sulted in invariance of the structure, the large density 
changes applied here (a factor of 10 and lO'^ respectively) 
do lead to systematic changes in the structure along the 
proposed isomorph, especially in the region of the first 
peak of g{f). This is not unexpected, because the in¬ 
crease in the steepness of the potential with decreasing 
density (see Fig. 6 ) makes close encounters between par¬ 
ticles less probable. This results in a steeper initial slope 
of g{f) and thus a higher first peak if the total number of 
nearest neighbors is unchanged [73] . Taking into account 
that the direct-isomorph-check isomorph “only” involves 
densities variations covering a single decade, we note that 
the deviation from isomorphic invariance is smallest for 
the state points obtained using A = 1. 

Not only the height of the peak, but also its position 
changes somewhat along the prospective isomorphs. This 
is an effect of the large change in the effective steepness 
of the potential (i.e., relative to ksT). The Yukawa po¬ 
tential has for instance both a BCC and an FCC solid 
phase with the triple point at k = 6.90 [19]. These two 
crystal structures have different reduced nearest-neighbor 
distances, and it is not surprising that some of these 


structural differences continue to exist in the liquid state. 
The change of the nearest-neighbor distance with density, 
however, indicates that a single value for A is not enough 
for accurately tracing out the isomorph. 


VI. DISCUSSION 

We have shown that the Yukawa system is Roskilde 
simple. This was shown from the strong correlations 
between equilibrium virial and potential-energy fluctu¬ 
ations and by the fact that the Yukawa system has iso¬ 
morphs in its phase diagram. An isomorph is a curve 
of constant excess entropy, but the curve can also be 
estimated directly from the pair potential. We have ver¬ 
ified that the dynamics are invariant on the isomorphs 
as predicted by the theory. The structure of the fluid 
as characterized by the radial distribution function was 
found to be invariant to a lesser degree, in part because 
the position of the first peak shifts slightly. Our estimate 
of the isomorph shape from the pair potential uses the 
nearest-neighbor distance. Consequently, Eq. (27) with a 
constant A can not give an exact description of the den¬ 
sity dependence of 7 along an isomorph. In view of the 
invariant dynamics shown in Fig. 9(a) and (b), however, 
Eq. (25) with A = 1 must be said to give a rather good 
description of the isomorph shape. 

In summary, we have shown that the isomorph the¬ 
ory provides a simple and general framework for under¬ 
standing previous findings for the Yukawa system. Thus 
our results explain the previously identified expression 
of the melting line of Yukawa systems [20, 21], which 
here follows from the fact that the melting line is an iso¬ 
morph [54]. Likewise, the isomorph theory [12, 28, 29] 
explains the observation that the transport coefficients 
of the Yukawa fluid collapse when plotted versus tem¬ 
perature scaled by the melting temperature (compare 
Eq. (16)), as well as the recently proposed expression 
for the ideal glass transition line of mode coupling the¬ 
ory [18] .Finally, we point out that the Yukawa system 
belongs to the “exponential” class of quasiuniversal sys¬ 
tems [62]. 


ACKNOWLEDGMENTS 

The centre for viscous liquid dynamics “Glass and 
Time” is sponsored by the Danish National Research 
Foundation via Grant No. DNRF61. 


[1] J. Rowlinson, Physica A 156 , 15 (1989). [3] P. Debye and E. Hiickel, Phys. Zeitschrift 24 , 185 (1923). 

[2] H. Yukawa, Proc. Physico-Mathematical Soc. Japan 17 , [4] B. Derjaguin and L. Landau, Acta Physicochim. U.R.S.S. 

48 (1935). 14 , 633 (1941). 














11 


[5] E. J. W. Verwey and J. T. G. Overbeek, Theory of the 
stability of lyophobic colloids (Elsevier, 1948). 

[61 A.-P. Hynninen and M. Dijkstra, Phys. Rev. E 68 , 
021407 (2003). 

[7] S. A. Khrapak, A. V. Ivlev, and G. E. Morfill, Phys. 
Plasmas 17, 042107 (2010). 

[ 8 ] V. E. Fortov, A. V. Ivlev, S. A. Khrapak, A. G. Khrapak, 
and G. E. Morfill, Phys. Rep. 421, 1 (2005). 

[9] G. E. Morfill and A. V. Ivlev, Rev. Mod. Phys. 81, 1353 
(2009). 

[10] P. N. Pusey and W. van Megen, Nature 320, 340 (1986). 

[11] H. Ikezi, Phys. Fluids 29, 1764 (1986). 

[12] M. O. Robbins, K. Kremer, and G. S. Grest, J. Chem. 
Phys. 88 , 3286 (1988). 

[13] E. J. Meijer and D. Frenkel, J. Ghem. Phys. 94, 2269 
(1991). 

[14] R. T. Farouki and S. Hamaguchi, Appl. Phys. Lett. 61, 
2973 (1992). 

[15] J. H. Ghu and L. I, Phys. Rev. Lett. 72, 4009 (1994). 

[16] H. M. Thomas, G. E. Morfill, V. Demmel, J. Goree, 
B. Feuerbacher, and D. Mohlmann, Phys. Rev. Lett. 73, 
652 (1994). 

[17] K. N. Dzhumagulova, T. S. Ramazanov, and R. U. 
Masheeva, Phys. Plasmas 20, 113702 (2013). 

[18] A. Yazdi, A. V. Ivlev, S. A. Khrapak, H. M. Thomas, 
G. E. Morfill, H. Lowen, A. Wysocki, and M. Sperl, Phys. 
Rev. E 89, 063105 (2014). 

[19] S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin, Phys. 
Rev. E 56, 4671 (1997). 

[20] O. S. Vaulina and S. A. Khrapak, J. Exp. Theor. Phys. 
90, 287 (2000). 

[21] O. S. Vaulina, S. A. Khrapak, and G. E. Morfill, Phys. 
Rev. E 66 , 016404 (2002). 

[22] F. A. Lindemann, Phys. Zeitschrift 11, 609 (1910). 

[23] J. J. Gilvarry, Phys. Rev. 102, 308 (1956). 

[24] T. S. Ingebrigtsen, T. B. Schrpder, and J. G. Dyre, Phys. 
Rev. X 2, 011011 (2012). 

[25] J. C. Dyre, J. Phys. Chem. B 118, 10007 (2014). 

[26] K. Kremer, G. S. Grest, and M. O. Robbins, J. Phys. A. 
Math. Gen. 20, L181 (1987). 

[27] H. Lowen and G. Szamel, J. Phys. Condens. Matter 5, 
2295 (1993). 

[28] Y. Rosenfeld, Phys. Rev. E 62, 7524 (2000). 

[29] Y. Rosenfeld, J. Phys. Condens. Matter 13, L39 (2001). 

[30] H. Ohta and S. Hamaguchi, Phys. Plasmas 7, 4506 

( 2000 ). 

[31] J. Daligault, Phys. Rev. E 86 , 047401 (2012). 

[32] S. A. Khrapak, O. S. Vaulina, and G. E. Morfill, Phys. 
Plasmas 19, 034503 (2012). 

[33] G. Faussurier, S. B. Libby, and P. L. Silvestrelli, High 
Energy Density Phys. 12, 21 (2014). 

[34] U. R. Pedersen, T. Christensen, T. B. Schrpder, and J. C. 
Dyre, Phys. Rev. E 77, 011201 (2008). 

[35] U. R. Pedersen, N. P. Bailey, T. B. Schrpder, and J. C. 
Dyre, Phys. Rev. Lett. 100, 015701 (2008). 

[36] N. P. Bailey, U. R. Pedersen, N. Gnan, T. B. Schrpder, 
and J. C. Dyre, J. Chem. Phys. 129, 184507 (2008). 

[37] A. A. Veldhorst, L. Bphling, J. C. Dyre, and T. B. 
Schrpder, Eur. Phys. J. B 85 (2012). 

[38] A. K. Bacher and J. C. Dyre, Colloid Polym. Sci. 292, 
1971 (2014). 

[39] T. S. Ingebrigtsen, T. B. Schrpder, and J. C. Dyre, J. 
Phys. Chem. B 116, 1018 (2012). 

[40] A. A. Veldhorst, J. C. Dyre, and T. B. Schrpder, J. Chem. 


Phys. 141, 054904 (2014). 

[41] L. Separdar, N. P. Bailey, T. B. Schrpder, S. Davatol- 
hagh, and J. C. Dyre, J. Chem. Phys. 138, 154505 (2013). 

[42] D. E. Albrechtsen, A. E. Olsen, U. R. Pedersen, T. B. 
Schrpder, and J. C. Dyre, Phys. Rev. B 90, 094106 
(2014). 

[43] A. Malins, J. Eggers, and C. P. Royall, J. Chem. Phys. 
139, 234505 (2013). 

[44] E. H. Abramson, J. Phys. Chem. B 118, 11792 (2014). 

[45] A. Henao, S. Pothoczki, M. Canales, E. Guardia, and 

L. C. Pardo, J. Mol. Liq. 190, 121 (2014). 

[46] S. Pieprzyk, D. M. Heyes, and a. C. Braka, Phys. Rev. E 
90, 012106 (2014). 

[47] S. Prasad and C. Chakravarty, J. Chem. Phys. 140, 
164501 (2014). 

[48] J. W. P. Schmelzer and T. V. Tropin, J. Non. Cryst. 
Solids 407, 170 (2015). 

[49] U. Buchenau, J. Non. Cryst. Solids 407, 179 (2015). 

[50] D. M. Heyes, D. Dini, and A. C. Braka, Phys. Status 
Solidi B (2015). 

[51] T. S. Ingebrigtsen, J. R. Errington, T. M. Truskett, and 
J. C. Dyre, Phys. Rev. Lett. Ill, 235901 (2013). 

[52] T. S. Ingebrigtsen and J. C. Dyre, Soft Matter 10, 4324 
(2014). 

[53] E. Lerner, N. P. Bailey, and J. C. Dyre, Phys. Rev. E 90, 
052304 (2014). 

[54] N. Gnan, T. B. Schrpder, U. R. Pedersen, N. P. Bailey, 
and J. C. Dyre, J. Chem. Phys. 131, 234504 (2009). 

[55] C. M. Roland, S. Hensel-Bielowka, M. Paluch, and 
R. Casalini, Reports Prog. Phys. 68 , 1405 (2005). 

[56] L. Bphling, T. S. Ingebrigtsen, A. Grzybowski, 

M. Paluch, J. C. Dyre, and T. B. Schrpder, New J. Phys. 
14, 113035 (2012). 

[57] Y. Rosenfeld, Phys. Rev. A 15, 2545 (1977). 

[58] T. S. Ingebrigtsen, A. A. Veldhorst, T. B. Schrpder, and 
J. C. Dyre, J. Chem. Phys. 139, 171101 (2013). 

[59] Y. Rosenfeld and P. Tarazona, Mol. Phys. 2, 141 (1998). 

[60] K. Y. Sanbonmatsu and M. S. Murillo, Phys. Rev. Lett. 
86 , 1215 (2001). 

[61] J. Daligault, Phys. Rev. Lett. 96, 065003 (2006). 

[62] A. K. Bacher, T. B. Schrpder, and J. C. Dyre, Nat. Com- 
mun. 5, 5424 (2014). 

[63] T. B. Schrpder and J. C. Dyre, J. Chem. Phys. 141, 
204502 (2014). 

[64] M. P. Allen and D. J. Tildesley, Computer Simulations 
of Liquids (Clarendon, Oxford, 1987), 1st ed. 

[65] J. C. Dyre, Phys. Rev. E 88 , 042139 (2013). 

[ 66 ] T. S. Ingebrigtsen, L. Bphling, T. B. Schrpder, and J. C. 
Dyre, J. Chem. Phys. 136, 061102 (2012). 

[67] N. P. Bailey, U. R. Pedersen, N. Gnan, T. B. Schrpder, 
and J. C. Dyre, J. Chem. Phys. 129, 184508 (2008). 

[68] Roskilde University Molecular Dynamics, 
http://rumd.org. 

[69] S. Toxvaerd and J. C. Dyre, J. Chem. Phys. 134, 081102 

( 2011 ). 

[70] J. P. Hansen, Phys. Rev. A 8 , 3096 (1973). 

[71] K.-C. Ng, J. Chem. Phys. 61, 2680 (1974). 

[72] W. L. Slattery, G. D. Doolen, and H. E. DeWitt, Phys. 
Rev. A 21, 2087 (1980). 

[73] L. Bphling, N. P. Bailey, T. B. Schrpder, and J. C. Dyre, 
J. Chem. Phys. 140, 124510 (2014). 

[74] N. P. Bailey, L. Bqhling, A. A. Veldhorst, T. B. Schrqder, 
and J. C. Dyre J. Chem. Phys. 139, 184506 (2013). 

[75] S. A. Khrapak and H. M. Thomas, Phys. Rev. E 91, 


12 


023108 (2015). [77] S. A. Khrapak, Phys. Plasmas 20, 054501 (2013). 

[76] S. A. Khrapak, N. P. Kryuchkov, S. O. Yurchenko, and [78] S. M. Stishov, JETP Lett. 67, 90 (1998). 

H. M. Thomas, J. Chem. Phys. 142, 194903 (2015). 


